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ABSTRACT 

Feedback by active galactic nuclei (AGN) appears to be critical in balancing radiative cooling of the low-entropy 
gas at the centers of galaxy clusters and in mitigating the star formation of elliptical galaxies. New observations of 
M87 enable us to put forward a comprehensive model for the physical heating mechanism. Low-frequency radio 
observations by LOFAR revealed the absence of fossil cosmic ray (CR) electrons in the radio halo surrounding 
M87. This puzzle can be resolved by accounting for the CR release from the radio cocoons and the subsequent 
mixing of CRs with the dense ambient intracluster gas, which thermalizes the electrons on a timescale similar to 
the radio halo age of 40 Myrs. Hadronic interactions of similarly injected CR protons with the ambient gas should 
produce an observable gamma-ray signal in accordance with the steady emission of the low state of M87 detected 
by Fermi and H.E.S.S. Hence, we normalize the CR population to the gamma-ray emission, which shows the same 
spectral slope as the CR injection spectrum probed by LOFAR, thereby supporting a common origin. We show 
that CRs, which stream at the Alfven velocity with respect to the plasma rest frame, heat the surrounding thermal 
plasma at a rate that balances that of radiative cooling on average at each radius — thereby setting the thermal 
pressure profile equal to that of the streaming CRs. However, the resulting global thermal equilibrium is locally 
unstable and allows for the formation of the observed cooling multi-phase medium through thermal instability. 
Provided CR heating balances cooling during the emerging "cooling flow", the collapse of the majority of the gas 
is halted at 1 keV, which is in accordance with x-ray data. We show that both the existence of a temperature floor 
and the similar radial scaling of the heating and cooling rates are generic predictions of the CR heating model. 
Subject headings: cosmic rays — galaxies: active — galaxies: clusters: intracluster medium — galaxies: indi- 
vidual (M87) — gamma rays: galaxies: clusters — radio continuum: galaxies 
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1. INTRODUCTION 

The central cooling time of the x-ray emitting gas in galaxy 
clusters and groups is bimodally distributed; approximately 
half of all systems have radiative cooling times of less than 
1 Gyr and establish a population of cool core (CC) clusters 
dCavagnolo et al. 2009). In the absence of any heating pro- 
cess, these hot gaseous atmospheres are expected to cool and 
to form stars at rates u p to several hundred M yr" 1 (see 
iPeterson & Fabianl 120061 for a review). Instead, the observed 
gas cooling and star formation rates are reduced to levels 
substantially below those expected from unimpeded cooling 
flows. High-resolution Chandra and XMM-Newton observa- 
tions show that while the gas temperature is smoothly declining 
toward the center, there is a lack of emission lines from gas at a 
temperat ure below about one third of the ambient cluste r tem- 
perature (IPeterson et al.l 120031 IPeterson & Fabianl 120061) . Ap- 
parently, radiative cooling is offset by a yet to be identified 
heating process that is associated with the AGN jet-inflated ra- 
dio cocoons that are co-localized with the cavities seen in the x- 
ray maps. The interplay of cooling gas, subsequent star forma- 
tion and nuclear activity appears to be ti ghtly coupled to a self- 
regul a ted fee dback loop (for reviews, see lMcNamara & Nulsenl 
I2007L 120121) . Feedback is considered to be a crucial if not 
unavoidable process of the evolution of galaxies and super- 
massi ve black holes dCroton et alj|2006t iPuchwein & Springell 
120131) . While the energetics of AGN feedback is sufficient to 
balance radiative cooling, it is unable to transform CC into 
non-CC clusters on the buoyancy timescale due to the weak 
coupling between the mechanical energy to the cluster gas 
dPfrommer et al.l 120121) . Matching this observational fact has 
proven to be difficult in cosmological simulations of AGN 



feedb ack that include realistic metal-line cooling dDubois et alj 
l20Tl . 

Several physical processes have been proposed to be respon- 
sible for the heating, includi ng an inward transport of heat from 
the hot outer c luster regions (INaravan & Medvedevl2001l) . tur- 
bulent mixing dKim & Naravanl2003l) . red istribution of heat by 
buoyancy-induced turbul ent convection (IChandran & Raseral 
120071: ISharma et aD2009l) . and dissipation of mechanical heat- 
ing by outflo ws, cocoons, or sound waves from the cen- 
tral AGN (e.g..lChurazov et al.l | 200lHBruggen & Kaiserjl2002l: 
iRuszkowski & Begelmanll2002l: iRuszkowski et alj|2004l) . The 
total energy required to inflate a cavity of volume V in an 
atmosphere of thermal pressure P t h is equal to its enthalpy 
£cav = 7/(7 _ l)PthV - 4.fthV, assuming a relativistic equation 
of state with 7 = 4/3 within the cocoons. Hence, the dissi- 
pation of mechanical energy resulting from the volume work 
done on the surroundings only makes use of one quarter of the 
total available enthalpy. Only a fraction thereof is dissipated 
in the central regions, which have the smallest cooling time 
(since most of the temperature increase caused by weak shocks 
is provided adiabatically). 

It would be more economic to additionally use the remain- 
ing part that is stored in the internal energy of the (presum- 
ably) relativistic particle population. If those CRs were mixed 
into the ambient gas, there would be a promising process to 
transfer CR into thermal energy. Fast streaming CRs along 
the magnetic fie ld excite Alfven waves through the "stream- 
ing instability" dKulsrud & Pearcdll969l) . Scattering off this 
wave field effectively limits the bulk speed of CRs. Damp- 
ing of these waves transfers CR energy and momentum to the 
thermal gas at a rate that scales with the CR pressure gradient. 
Hence, this could provide an efficient means of suppressing 
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the cooling catastrophe of CCs JGuo & Ohll2008t lEnfilin et al] 
HoTltlFuiita & Ohirall2012HWiener et al. 1120 131) . but only if the 
CR pressure gradient is sufficiently large. 

This work studies the consequences of this CR-heating 
model with the aim of putting forward an observationally well- 
founded and comprehensive model. In particular, we intend 
to eliminate two of the most uncertain assumptions of this 
process, the mixing of CRs with the ambient plasma and the 
CR abundance in the central cluster regions. Moreover, we 
tackle the question about thermal balance in the core region 
of a CC cluster, but constrained to this model in which CR 
Alfven wave heating competes with radiative cooling. To this 
end, we focus on the closest active galaxy that is interacting 
with the cooling cluster gas that provides us with the most de- 
tailed view on AGN feedback possible. This selects M87, the 
cD galaxy of the Virgo cluster at a distance of 16.7 ± 0.2 Mpc 
dMei et al.ll2007h . where 1' corre sponds to 4.85 kpc. It is well- 
studied not only at x-rays (e.g.. iMillion et al] 12010) but also 
and especially at non-thermal wave length s, i.e., in the radio 
( Owen et"aT]|2000t Ide Gasperin et al.l 12012b and 7-ray regime 
dRieger & Aharonianl2012l) . 

In Section|2]we show how new observations of M87 by LO- 
FAR, Fermi, and H.E.S.S. can shed light on the injection and 
mixing of CRs into the ambient plasma and estimate the CR 
pressure contribution in the central regions. Section|3]discusses 
the balance of CR heating with radiative cooling. Particular at- 
tention is paid to local thermal stability of the global thermal 
equilibrium that we find. In Section [4] we present predictions 
of our model that include gamma-ray, radio, and Sunyaev- 
Zel'dovich observations and discuss how those can be used to 
quantify open aspects of the presented heating mechanism. Fi- 
nally, we compare our model with other solutions to the "cool- 
ing flow problem" and conclude in Section [5] In Appendix lAl 
we derive a thermal instability criterion for a fluid that is cool- 
ing by thermal bremsstrahlung and line emission and is heated 
by a generic mechanism while allowing for thermal and CR 
pressure support. 

2. A NON-THERMAL MULTI-MESSENGER APPROACH 

2.1. Cosmic-Ray Injection: Insights by LOFAR 

The characteristic frequency v s at which a CR electron (CRe) 
emits synchrotron radiation scales with the magnetic field 
strength B and the CRe's cooling time r as v s oc B~ 3 t~ 2 (assum- 
ing the strong-field limit for which inverse Compton cooling 
is negligible). Hence, observations at lower radio frequencies 
should enable to probe older fossil CRe populations. Those are 
expected to accumulate at lower energies, which would imply 
a low-frequency spectral steepening. Equivalently, such obser- 
vations should unveil the time-integrated non-thermal action of 
AGN feedback that is expected to be spread out over larger 
spatial scales. However, recent LOFAR observations of Virgo 
A, the radio source of M87, provided two unexpected findings. 
(1) The observations do not reveal a previously hidden spec- 
tral steepening toward low frequencies and show the extended 
radio halo to be well c onfined within boundari es that are identi- 
cal at all frequencies (Ide Gasperin et al.ll2012l) . The possibility 
that this may point to a special time, which is characterized by 
a recent AGN outburst after long silence, is unlikely because 
statistical studies of the radio-mode feedback in complete sam- 
ples of clusters imply that the duty cycle of AGN outbursts 
with the potential to heat the gas significantly in CC clusters is 
> 60% and could approach unity after accounting for selection 
effects dBirzan et alj|20l2l) . (2) While the core region displays 
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Figure 1. Cooling time of electrons for typical conditions in the central region 
of M87 as a function of normalized electron momentum, where p is the physi- 
cal momentum. During their advective transport in the radio cocoon, the (high- 
energy) electrons experience inverse Compton (IC), synchrotron, and adiabatic 
losses, for which we assume an expansion factor of / = 10. After the electrons 
mix with the ambient thermal plasma, Coulomb interactions can thermalize 
them. Electrons with kinetic energy E ~ pm c c~ ~ 150 MeV are the most 
long-lived, with a maximum cooling time of r ~ 40 Myr that corresponds to 
the radio halo age. 

a radio spectral indexQof a v = 0.6 ±0.02 (consistent with an in- 
jection CR spectral index of a = 2a u + \ = 2.2 ±0.04), the low- 
frequency halo spectrum outside the central c ocoon is consider- 
ably st eeper with a„ = 0.85. As suggested by de Gas perin et al.l 
d2012l) . the most probable explanation of this steepening at low 
frequencies is strong adiabatic expansion in combination with 
a superposition of differently synchrotron-aged CRe popula- 
tions. Adiabatic expansion across a density expansion factor 
/ = Pij P\ from the compact source to the more dilute halo re- 
gions shifts all CR momenta by the factor pi = /'' 3 /?L CRes 
that are leaving the cocoon regions have spent different time 
intervals in regions of strong magnetic fields, which leads to a 
distribution of spectra with different break momenta. Superpo- 
sition of those spectra steepens the spectrum down to low fre- 
quencies if the cooling time of the most strongly cooled CRes 
has been comparable to the transport time from the core to the 
current emitting region. 

In Fig. [D we sh ow the cooling time of CRes (see 
iPfrommer et al.l 120081 for definitions). CRes of energy ~ 
200 (400) MeV that are responsible for - 25 (100) MHz ra- 
dio emission in a 20yuG field have energy loss timescales of 
~ 150 (75) Myr, which drops further by the factor f 1 ' 3 after 
accounting for adiabatic expansion. CRes can cool through 
inverse Compton (IC) interactions with various seed photon 
fields, including the cosmic microwave background (CMB), 
star light, and dust-reprocessed infra-red emission. Within 
the central radio halo region (r < 34 kpc), the star formation- 
induced photon fiel ds at least contribute as much as the CMB 
dPinzke et al.l 1201 lh . which motivates our choice of a seed 
photon energy density that is twice that of the CMB. Syn- 
chrotron cooling in equipartition magnetic fields of 10-20/iG 

1 We define spectral indexes according to F v oc v~ a " . 
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(Ide Gasperin et al. 2012) and adiabatic expansion substantially 
decreases the cooling times further during the advective trans- 
port of CRes in the radio cocoons. Those rise in the cluster 
atmosphere (either buoyantly or d riven by the jet) an d uplift 
dense gas from the central regions (IWerner et al.lEoTol) . Those 
rising motions of light relativistic fluid cause strong downdrafts 
at the sides of the cocoons, which induce Rayleigh-Taylor and 
Kelvin-Helmholtz instabilities that shred the cocoons and ex- 
cite turbulence, which mixes the CRes with the dense uplifted 
ambient plasma. This picture appea r s to be supported by ra- 
dio maps of M87 (lOwen et al.ll2000t Ide Gasperin et al.ll2012l) 
that show a diffusely glowing radio halo surrounding the co- 
coons, which themselves show the characteristic mushroom 
shape of a Rayleigh-Taylor instability. Coulomb interactions of 
the released CRes with the dense ambient intracluster medium 
(ICM) cause them to be thermalized after a maximum cooling 
time of t ~ 40 Myr, which corresponds to the radio halo age 
(see Fig.[TJ. Thus, this naturally explains the absence of a fos- 
sil CRe signature in the low-frequency radio spectra as well as 
the spatially similar morphology of Virgo A at all frequencies. 

Diffusive shock acceleration operates identically on relativis- 
tic particles of the same rigidity (R = pc/ze, where p is the 
physical momentum and z is the ion charge) and unlike elec- 
trons, protons do not suffer an injection problem into the ac- 
celeration process and have negligible radiative losses. Hence 
we expect CR protons (CRps) to dominate the pressure of 
the cocoon£| and to exhibit a similar injection spectrum as 
the CRes in the central regions (a « 2.2). Similarly to the 
CRes, CRps should also be released from the cocoons and 
mixed into the ambient ICM. Since adiabatic expansion leaves 
their spectral shape invariant and CRps have negligible radia- 
tive losses, we do not expect a spectral steepening of the CRp 
population in the halo. A necessary consequence of this pic- 
ture are hadronic interactions of CRps with the ambient gas. 
This should produce a ste ady gamma-ray signal from M87 
dPfrommer & En Blin 2003), with a spectral index above GeV 
energies of a 7 = a rj 2.2. 

2.2. Cosmic-Ray Pressure Estimates from Fermi and H.E.S.S. 

To test this hypothesis, we turn to the gamma-ray regime. 
A joined unbinned likelihood fit with a power-law (dN/dE oc 
E~ a i) to the Fermi and H.E.S.S. data of M87's low state 
yields F(> 1 GeV ) = 1.35(±0 .35) x 10" 9 photons cm" 2 s" 1 and 
a 7 = 2.26 ±0.13 (lAbdo et al.ll2009l) . which is consistent with 
our expectation if the gamma-rays are produced in hadronic 
CRp-p interactions. Moreover, a comparison of the Fermi 
light curve for the last 14 months with the initial 10-months 
data r evealed no evidence of variability in the flu x above 100 
MeV dAbramowski et al.ll2012t lAbdo et al.ll2009l)FI Gamma- 
ray upper limits on M87 set by EGRET (iReimer et all 120031) 
are consistent with the steady flux seen by Fermi. This demon- 
strates that the EGRET and Fermi observations of M87, which 
together span more than two decades, show no apparent in- 
crease in flux. Additionally, the low flux state at very high en- 
ergies E > 100 GeV d oes not provide any hint of variability 
dAharonian et alj|2006l) . although the low signal statistics ren- 

2 In fact, a-particles cany a significant fraction of the total CR energy, which 
we absorb into the proton spectrum with the following line of reasoning. A 
GeV energy a-particle can be approximated as an ensemble of four individ- 
ual nucleons travelling together due to the relatively weak MeV nuclear bind- 
ing energies in com parison to the kinetic energy of relativistic protons. See 
EnBlin et al. 12007) for an extended discussion. 

3 A light curve analysis of 4 years of Fermi data does not show any sign of 
variability either (F. Rieger, priv. comm.). 



ders such a test very difficult. Encouraged by the spectral sim- 
ilarity of radio and gamma-ray observations and the constant 
gamma-ray flux, we will now assume that the gamma rays are 
tracing hadronic CR interactions. 

To estimate the CR pressure that is implied by the ob- 
served gamma-ray emission, we choose a constant CR-to- 
thermal pressure ratio, Xqr = Pcr/Poi (which we require for 
self-consistency as we will see in Section[3]). We assume a CR 
power-law population, 



d 4 N 
dpd 3 x 



■ Cp- a 6(p-q), 



(1) 



where 9 denotes the Heaviside step function and a = 0.226. We 
decided in favor of the spectral index of the gamma-ray emis- 
sion, which likely probes the prevailing CRp distribution in the 
radio halo region. In contrast, the radio emission of the central 
cocoon region may signal current acceleration conditions that, 
albeit generally very similar, may be slightly different from the 
ones that have accelerated the CRs in the past, which have now 
been released. 

We use a fit to the spherically symmetric thermal pressure 
profile as inferred from x-ra y observations ([Matsushita et al.l 
12001 IBohringer et al.l [T994b . The thermal pressure is Pth = 
nkT = n e kT / ' p, e , where fj, e = n e Jn = fi(Xu + l)/2 ~ 0.517 for 
a fully ionized gas of primordial element composition with a 
hydrogen mass fraction X H = 0.76 and mean molecular weight 
\i = 0.588. The electron density and temperature profilesQof the 
central region of the Virgo cluster (r £ 100 kpc corresponding 
to 20') are given by fits to the x-ray profiles and have the form 



n e (r) = 



T e (r) = T + (Ti-T ) 



-3/3/ 



1/2 



and 



1 + 



(2) 
(3) 



where (n u n 2 ) = (0.15,0.013) cm" 3 , (n,r 2 ) = (1.6,20) kpc, 
Whfc) = (0.42,0.47), (kT ,kTi) = (1,3) keV, and r, = 13.5 kpc. 

The sharp drop of radio morphology of Virgo A's halo sug- 
gests that the CRps are also confined to a central region of 
radius r max = 34 kpc. Hence we integrate the gamma-ray 
flux due to hadronic CRp-p interactions from within a spher- 
ical region bounded by r^ax- W e use a semi-analytic formal- 
ism ( Pfr ommer & E nBlin 2004) that ha s been modified to con - 
strain the CR-to-thermal pressure ratio dPfrommer et al.ll2008l) . 
We obtain Xcr = 0.38 for a vanishing low-momentum cutoff, 
q = 0. Accounting for a non-zero low-momentum cutoff of 
the CRp distribution as a result of Coulomb cooling, we ob- 
tain X C r = 0.31 (0.34) for q = 0.58 (0.25) which corresponds 
to a Coulomb cooling time of CRps of r = 40 Myr, assuming 
an ambient electron density of n e = 10 _1 (10~ 2 ) cm" 3 . Since 
there is considerable uncertainty about the shape of the low- 
momentum regime of the CR distribution function, we decided 
to adopt the most conservative value of Xqr = 0.31 for the re- 
minder of this work. We note that a fresh supply of CRps could 
lower q and therefore increase the CR pressure. 

We checked that the radio emission from secondary CRes, 
which are inevitably injected by hadronic CRp-p reactions, has 

4 In Virgo, the gas is observed to have a maximum temperature T ~ 
3 X 10 7 K at a clust er-centric radius of ~ 200 kpc (corresponding to 40', 
IBohringer etal.119941) . 



4 



PFROMMER, C. 



a subdominant contribution to the observed radio emission. 
While the total radio emission due to secondaries is well be- 
low the observed flux density of the radio halo of M87, fu- 
ture high-frequency observations at £ 10 GHz of the external 
halo parts may be abl e to p robe the hadronically-induced radio 
emission (see Section 14721 for more details). At these frequen- 
cies, the hadronically-induced emission component cannot any 
more be neglected and complicates inferences on the CRe pitch 
angle distribution in different spectral ageing models of radio- 
emitting CRes. 

In general, cluster merger and accretion shocks are also ex- 
pected to accelerate CRs, which add to the CR population in- 
jected by the central AGN. Non-observations of gamma-ray 
emission at GeV and TeV energies limits their contribution 
to Xq r < 0.017 in the Com a and Perseus clusters (lArlen et al.l 
120121 ; lAleksic et alj 120 12h and typic ally to less than a few 
per cent for the n ext to best targets (lAckermann et al.l 120101 
iPinzke et al.ll2.Ql lb . We hence neglect their contribution in the 
center, where the AGN-injected population apparently domi- 
nates. 

These findings for Xcr are in agreement with constraints on 
the non-thermal pressure contribution in M87 that have been 
derived by comparing the gravitational potentials inferred from 
stellar and globular cluster kinematics and fro m assuming hy- 
drost atic equilibrium of the x-ray emitting gas (IChurazov et alJ 
120101) . Depending on the adopted estimator of the optical 
velocity dispersion, they find a non-thermal pressure bias of 
X at = P M /P th « 0.21-0.29 (which is even higher in M87 for 
more sophisticated estimates that are less sensitive to the un- 
known orbit anisotropy). There are several effects that could 
(moderately) increase X nt . (1) A more extended distribution 
of P nt in comparison to P t h could hide a larger non-thermal 
pressure contribution (since only the pressure gradient enters 
into the equation for hydrostatic equilibrium), (2) the pres- 
ence of an x-ray "invisible" cool gas phase that is coup led to 
the hot phase could bias X nt low (IChurazov et afl 120081) . and 
(3) the tension in the orbit anisotropy estimates of the stars 
and globular clusters may call f or triaxial potential models 
dRomanowskv & Kochanekll 20011) . 

2.3. Pressure Content of Radio Cocoons 

So far, we implicitly assumed that the pressure of radio 
cocoons is provided by a relativistic component. While it 
is unknown what physical component provides the dominant 
energetic contribution to radio cocoons, here we will review 
the reasons for our assumption. Equipartition arguments for 
the radio emitting cocoons show that the sum of CRes and 
magnetic fields can only account for a pressure fraction of 
~ 10% in comparison to the external pressure, with which 
the coco ons are apparently i n approximate hydrosta tic equi- 
librium (Bla nton et al.l 12003; de Gasperin e t al. 2012). How- 
ever, these estimates rely on power-law extrapolations from 
the high-energy part of the CRe distribution (with Lorentz fac- 
tor 7 ^ 10 3 , which was visible so far at radio frequencies 
v >, 300 MHz) to lower energies which dominate the CRe en- 
ergy budget, leaving the possibility of an aged fossil CRe popu- 
lation that would dominate the pressure. LOFAR observations 
have now unambiguously demonstrated the absence of such a 
hypot hetical fossil CRe population in M87 (de Gas perin et al.l 
120121) . Since CRps should be accelerated at least as efficiently 
as the visible CRes and since CRps have negligible radiative 
cooling rates (in comparison to CRes), we make the very plau- 
sible assumption in this work that CRps provide the dominating 
pressure component of cocoons (and work out a few smoking- 



gun tests for this scenario in Section |4}. The alternative sce- 
nario of a hot, but subrelativistic (Maxwellian) proton popu- 
lation is unlikely because of three reasons. (1) CRp Coulom b 
cooling time scales are too short (see, e.g., EnBlin et al. 2007), 
implying a burst-like heating once the content has been set 
free as a result of shredding the radio cocoons by the action 
of Kelvin-Helmholtz and Rayleigh-Taylor instabilities. How- 
ever, there is no signature in the x-ray observati ons of such a 
locali zed heating around the shredded cocoons dMillion et al.l 
120101) . (2) Moreover, LOFAR observations show no sign of 
such a ^ 100 MeV electron distribution, requiring a helicity 
fine-tuning for the acceleration process. (3) Adiabatic expan- 
sion of a composite of relativistic and non-relativistic (thermal) 
populations favors the CR over the non-relativistic population 
since Pcr/Aii oc / _1 ^ 3 , due to the softer equation of state of 
CRs. Here we assumed an ultra-relativistic CR population with 
an adiabatic exponent of 7 = 4/3 and / = pz/ pi is the density 
expansion factor starting from the densities at which the two 
populations got energized. 

3. COSMIC-RAY HEATING VERSUS RADIATIVE COOLING 

The previous arguments provide strong evidence that CRs 
are released from the radio cocoons of M87 and mixed into the 
ICM on a time scale of r ~ 40 Myr. Interpreting the observed 
gamma-ray emission as a signature of hadronic CRp-p inter- 
actions enabled us to estimate the CR pressure contribution. 
We now turn to the question whether the resulting CR Alfven 
wave heating is powerful enough to balance radiative cooling 
and what this may imply for thermal stability of the gas in M87 
and CCs in general. 

3.1. Cosmic-Ray Heating 

Advection of CRs by the turbulent gas motions produces 
centrally enhanced profiles. In contrast, CR streaming along 
magnetic field lines implies a net outward flu x of CRs that 
causes the profiles to flatten (E nBlin et al.ll201 ll) . A necessary 
requirement for CR streaming to work is a centrally peaked 
profile of the CR number density. Do we ha ve any evi- 
dence for this in M87? Radi o maps of Virgo A dOwen et al.l 
2000; de Gas perin et al.l2012h suggest that turbulent advection, 
driven by the rising cocoons, is strong enough to maintain a 
sufficiently strong inward transport of CRs. Alternatively, the 
smoothly glowing radio halo that fills the central region for 
f < '"max may be caused by escaping CRs from the cocoons al- 
ready at the early stages of their rise. Following observational 
and theoretical arguments, the current distribution of the CR 
pressure appears to follow a centrally peaked profile. 

The CR Alfven wa ve heating term is given by dWentzell 
rTWTtlGuo^Orill2008h 



■Hcr = -v A VPcr = -va V r (P C r) n + 



SPt 



CR 



61 



(4) 



where v A = B/y/Aitp is the Alfven velocity, p is the gas mass 
density, and we decompose the local CR pressure gradient 
into the radial gradient of the azimuthally averaged CR pres- 
sure profile, (Pcr)u, and a fluctuating term, SPqr/SI = VPcr - 
V r (PcR)fi, which in general also points in non-radial direction. 
In steady state, CR streaming transfers an amount of energy to 
the gas per unit volume of 



Ae th = -T A V A VP C R « PCR = XcR^th, 



(5) 



where we used the Alfven crossing time t a = SI/va over the 
CR pressure gradient length, SI. Comparing the first and last 
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term of this equation suggests that a constant CR-to-thermal 
pressure ratio is a necessary condition if CR streaming is the 
dominant heating process, i.e., the thermal pressure profile ad- 
justs to that of the streaming CRs in equilibrium so that we will 
take the radial profile of P t h as a proxy for Pqr. 

The observed radial profile of P t h averages over the local 
pressure inhomogeneities within a given radial shell, especially 
toward the center where the azimuthally averaged radial pro- 
file plateaus out. However, such a constant P t n at the center 
is in strong disagree ment with the fluctu ating thermal pressure 
maps (see Fig. 7 of lMillion et alj|2010h . Those pressure fluc- 
tuations could have their origin in shocks, turbulence, pertur- 
bations of the gravitational potential, or entropy variations of 
advected gas that has not had time to come into pressure equi- 
librium with the surroundings. To illustrate one of these possi- 
bilities, we consider pressure fluctuations due to weak shocks 
and expand the Rankine-Hugoniot jump conditions in the limit 
of small Mach number, JA\ — 1 = s -C 1 : 



Pi 
P~i 



2 7 A^-(7-l) 
7+1 



1 + 



4 7 
7+ 1 



1 + ^, 



(6) 



where subscripts 1 and 2 denote upstream and downstream 
quantities in the shock rest frame and we assume 7 = 5/3 in 
the third step. Adopting a constant CR-to-thermal pressure ra- 
tio, we obtain an estimate for the CR heating due to pressure 
fluctuations as a result of a sequence of weak shocks, i.e., the 
second term on the right-hand side of Equation (|4) 
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Here we assume that the characteristic size scale of the pres- 
sure fluctuations is of order the cluster-centric radius, which 
can be realized, e.g., for spherically symmetric shocks whose 
curvature radius increases as r. We note, that some of the other 
mechanisms that may be responsible for the observed pressure 
fluctuations can in principle have a different scaling with radius 
and hence, have a modified radial dependence of the fluctuat- 
ing term in Hcr, in particular at the center. It will be subject to 
future work to quantify the contribution of those processes to 
the CR heating rate. 

We choose conservative values for the magnetic field 
strengths that are about half as much as t h e min imum energy 
estimates presented in Ide Gasperin et al.l (120121) . We adopt 
B = 20 (iG sfn e /0.1cm" 3 , which yields B = 24 (5) fiG for the 
central value (at the halo boundary, r max = 34 kpc). The density 
scaling yields a constant Alfven velocity v& ~ 130 km s -1 . The 
plasma beta factor j3 = P±/Pb increases from 1 1 (at the center) 
to 28 (at r max ), which corresponds to a magnetic pressure con- 
tribution of X B = 0.08 (center) and 0.035 (r max ). 

3.2. Global Thermal Equilibrium 

Radiative cooling of cluster gas is dominated by thermal 
bremsstrahlung and collisional excitation of heavy ions to 
metastable states and the subsequent line emission following 
de-excitation (the latter of which dominates for particle ener- 
gies 3 keV). These processes are characterized by the cool- 



ing rate 



C r ad = rc c «,A coo i(r,Z), 



(8) 



where the target density n, is the sum of ion densities and 
we defined the cooling function A coo i(r,Z) that depends on 
temperature T and metal abundance Z. The cooling func- 
tion is fairly flat between the range of temperatures of kTo = 
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Figure 2. CR heating versus radiative cooling rates. The cooling rates (blue) 
account for a fluctuating metal distribution in the range (0.7- 1.3) Zq and are 
globally balanced at each radius by CR Alfven wave heating. We show the 
CR heating rate Kcr due to the azimuthally averaged pressure profile that 
"smoothes" out the CR pressure gradient (dashed red) and the Wcr profile, 
where we additionally account for pressure fluctuations due to weak shocks of 
typical Mach number A4 = 1.1 (solid red). The CR heating rates become very 
uncertain outside the boundary of the radio halo of Virgo A (indicated with a 
dotted red). 



(1-2) keV, we encounter in the central regions, bounded by 
the radio halo. The X-ray inferred metallicity map shows sub- 
stantial spatial variations on small scales ranging in between 
Zq = (0.7 — 1.3) Zq, which are of similar a mplitude as the over- 
all radial trend in Z dMillion et al.l 120 10l> . Hence, we decided 
to bracket the cooling function for these ranges of tempera- 
ture and metallicity, which yields A coo i(7b,Zo) ps (1.8-2.8) x 
10~ 23 ergcm 3 s" 1 wh ile assuming collisional ion ization equi- 
librium (see Fig. 8 of Sutherl and & Dopitalll993l) . 

Figure|2]shows the radial profile of the radiative cooling rate 
C la d and that of the CR heating rate Hcr- For radii r £ 2 kpc, 
CR heating balances cooling already for the "smoothed" gra- 
dient pressure term (V, (Pcr)si)- The central pressure plateau 
causes the CR pressure gradient and hence, the CR heating 
to drop. To model the large visible fluctuations in P t h (and 
hence in Pqr according to Equation|5]l, we assume that they are 
caused by weak shocks of typical Mach number M. = 1 . 1 with 
a size scale that is of order the radius (i.e., we adopt e = 0.1 
in Equation [6}. This fluctuating pressure term appears to be 
critical in balancing the cooling rate with CR heating in the 
central region for r < 2 kpc. We conclude that CR Alfven 
wave heating can globally balance cooling, hence justifying 
our assumption of a constant CR-to-thermal pressure. Partic- 
ularly notable is the non-trivial result that the radial shapes of 
C ra d oc « 2 and 7-^cr oc V r P t h closely track each other. In fact, if 
CR heating dominates, P tri should be determined by the pres- 
sure profile of the streaming CRs and hence this resemblance of 
the radial heating and cooling rate profiles is an important sig- 
nature of self-consistency. The sharp boundaries of the radio 
halo of Virgo A may suggest that CRs are confined to within 
its boundaries, perhaps as a result of azimuthally wrapped mag- 
netic fields. Hence the CR heating rates become very uncertain 
outside the maximum halo radius. 

We will argue that this similarity of the CR heating and ra- 
diative cooling profiles is not a coincidence in M87, but rather 
generically expected in models of CR Alfven wave heating. 
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For simplicity, we use a constant gas-to-dark matter ratio in- 
side the dark-matter scale radius, i.e., n oc r , which is char- 
acteristic for the central density profile in CCs. Outside the 
temperature floor, we adopt a temperature profile with a small 
positive slope, T oc r a ' (with a, 0.3 typically) so that the 
pressure profile scales as P m = nkT oc r a '~ l . As a result, we ob- 
tain a scaling of the cooling and CR heating rates with radius of 
Crad r ~ 2 an( l ^cr oc d(r a '~ l ) / dr oc r a '~ 2 . Those are identical 
for a constant temperature (a, = 0). For a smoothly rising tem- 
perature profile (a, £ 0.3) the heating rate profile is tilted and 
slightly favoured over the cooling rates at larger radii, provided 
they match up at some inner radius (see Fig. |2). A corollary 
of this consideration is that the onset of cooling is smoothly 
modulated from the outside in, as we will see in the following. 



S 




'ocean of instability" 



3.3. Local Stability Analysis 

CR heating appears to be able to keep the gas in global ther- 
mal equilibrium in M87. However, this does not necessar- 
ily imply local thermal stability. We consider small isobaric 
perturbations to this equilibrium configuration, which main- 
tain pressure equilibrium at each radius. The assumption of 
isobaric perturbations is well justified for small timescales in 
comparison to the dynamical timescale, Td yn ps (2nGp m )~ x l 2 w 

150Myr x n\ , where «_2 = 10~ 2 cm~ 3 , p m = p b /f b is the total 
matter density, p b and f b w 0. 166 are the baryonic matter den- 
sity and baryon-to-dark matter mass fraction, respectively. In 
AppendixlAlwe derive a general instability criterion for a fluid 
that is cooling by thermal bremsstrahlung and line emission 
and heated by a generic mechanism with a power-law depen- 
dence on p and T while allowing for thermal and CR pressure 
support. We will start the discussion with a transparent order 
of magnitude calculation that provides insight into the physics 
and complements the formal treat ment in Appendi x lAl 

The generalized Field criterion (lFieldll9 65; Balb usll986l) for 
thermal instability in the presence of isobaric perturbations is 
given by 

d(C/T) 



10 s 



dT 



<0, 



(9) 



where C is defined as the net loss function per unit mass such 
that pC = C — H, where C and H are the cooling and heating 
rates per unit volume. On short timescales that we consider 
for our stability analysis, the CR population can be considered 
adiabatically since it does not suffer non-adiabati c losses such 
as hadronic interactions (e.g.. lEnfilin et al.l 12007b . In particu- 
lar, CR streaming i s an adiabatic process for the CR population 
(EnB lin et al.ll201 ll) . implying locally Pgr oc p 7 , where 7 is the 
adiabatic exponent of the CR population. We emphasize that 
this local adiabatic behavior of Pgr does not hold for the equi- 
librium distribution of CRs that is subject to non-adiabatic gain 
and loss processes, which modify the "CR entropy" (i.e., its 
probability density of the microstates)Q We obtain the scaling 
of the CR heating rate with temperature, 



Hcr = -v a VPcr oc p^ 3 oc T^~ l l\ 



(10) 



5 This is similar to thermal cluster gas that increases its pressure P t i, oc p 5 / 3 
upon adiabatic changes. However, non-adiabatic processes such as shocks trig- 
gered by gravitational infall, non-gravitational heating, and radiative cooling 
modify the gas entropy. As a result, the gas obeys an effective adiabatic ex- 
ponent 7 c ff = 91nP t h/91np 1.1 — 1.3 within a cluster that only slowly in- 
creases beyond the radii of accretion sh ocks toward the adiabatic value (e.g., 
IBattagli a et al. 2010; Capelo et al.120121) . 
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Figure 3. Local stability analysis of a fluid for which CR heating globally 
balances radiative cooling by bremsstrahlung and line emission (assuming so- 
lar metallicity). We show the instability criterion, T> oc d(C/T)/dT\ p , as a 
function of temperature for parameters appropriate for M87 (red) and for a 
ten times lower CR pressure (blue dashed). Precisely, we show arsinhCD), 
which has a linear (logarithmic) scaling for small (large) values of the argu- 
ment. Once the gas temperature drops below 3 X 10 7 K, it becomes thermally 
unstable. Assuming that the gas cools while maintaining approximate thermal 
equilibrium, its collapse is halted at around 10 7 K by the first "island of stabil- 
ity", which is consistent with the temperature floor seen in x-ray observations. 
However, the high-density tail of density fluctuations should able to cross this 
island and become again subject to thermal instability, likely sourcing some of 
the observable multi-phase gas. 

Here we identified the characteristic size of our unstable clump 
with an exponential pressure scale height h so that the lo- 
cal gradient length scales as 81 ~ h oc p~ 1/3 and we assumed 
Va = const. We suppress the dependence on pressure, which is 
constant for isobaric perturbations. The radiative cooling rate 
scales with temperature as 



c ad = P 2 \k z 2 t x i 2 +k,(t,z)\ oc t-^ 2 + v t- 2+x 



(11) 



where the first and second term in the brackets denote the 
bremsstrahlung and line-emission cooling. Here we introduce 
the logarithmic temperature slope of the line cooling function 
X and the ratio of line-to-bremsstrahlung cooling rate 77. For an 
ultra-relativistic CR population with 7 = 4/3, we find a stable 
fixed point in the bremsstrahlung regime since the slope of the 
heating rate (= -5/3) is steeper than that of the cooling rate 
(= -3/2). Increasing (decreasing) the temperature of a per- 
turbed parcel of gas implies a stronger cooling (heating) that 
brings the parcel back to the background temperature. How- 
ever in the line-cooling regime (at (1 -3) x 10 7 K), which is 
characterized by 77 ^ 1 and x ~ 0, the slope of the cooling rate 
( ,$ -2) is steeper than that of the CR heating rate (= -5/3), 
implying instability. 

This simple order of magnitude calculation is substantiated 
in Fi g. [3] w hich shows the instability criterion derived in Equa- 
tion (1A1 II) and restricted to CR Alfven wave heating as the sole 
source of heat in the central region of the Virgo cluster. A gas 
in thermal equilibrium ( pC = 0) is thermally unstable if 



(5/3-7X1+7?) 
1+X C r7 



<0. 



(12) 



We compute the logarithmic temperature slope of the line cool- 
ing function x(T) and the ratio of line-to-bremsstrahlung cool- 
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ing rate rj(T) for a gas at sol ar metallicity and assume c olli- 
sional ionization equilibrium (ISutherland & Dop ita 1993!). In 
Fig. [3] we show V as a function of temperature and adopt val- 
ues ap propr iate for M87 (7 = 1.37 and Xqr = 0.31, see Ap- 
pendix |A2]i. For most of the temperature range, we find that 
T> < and the gas is subject to thermal instability. However, 
there are a few "islands of stability" where T> > 0. 

The location of such an "island of stability" corresponds to 
a sufficiently positive slope of the cooling function A coo i as a 
function of T (i.e., it is located at the low-temperature side of 
every line-cooling complex). Increasing the metallicity Z in- 
creases the cooling rate and sharpens the line complexes. This 
implies more stable "islands of stability", i.e., they protrude 
more above the "ocean of instability", but they also extend over 
a slightly smaller temperature range. However, increasing Z 
does not significantly change their localization in temperature 
(as long as the CR heating can compensate the increasing cool- 
ing rate that accompanies the elevated Z). 

Apart from parameters associated with radiative cooling 
(e.g., metallicity), T> also depends on cluster-immanent param- 
eters associated with the CR distribution, i.e., the CR adiabatic 
exponent 7 and the CR-to-thermal pressure ratio Xqr. Of those, 
only Xqr can in principle vary substantially among clusters. To 
address the robustness of the shape of the instability criterion T> 
with respect to variations in Xcr, in Fig. Owe additionally show 
T> for a ten times lower CR pressure (which would however not 
suffice to maintain global thermal equilibrium in M87). Most 
importantly, the shape of T> versus temperature barely changes, 
except that the second "island of stability" vanishes for low val- 
ues of Xcr. 

3.4. A Critical Length Scale for the Instability 

Before we can discuss the fate of the cooling gas, we will be 
considering the typical length scale of a region that becomes 
subject to thermal instability. We expect cooling gas parcels 
only to appear in systems whose dimension is greater than a 
critical length scale A cr ; t , below which CR Alfven wave heat- 
ing smoothes out temperature inhomogeneities. We will esti- 
mate A cr ; t by considering thermal balance for a cool parcel of 
radius r embedded in a medium with an equilibrium CR pres- 
sure Pes. = -XcR-fth- CR streaming transfers an amount of en- 
ergy to a given parcel at a rate ~ r 3 f s VA(Pcs./r), where f s is a 
suppression factor that depends on the magnetic topology con- 
necting the parcel with its surroundings. Line emission and 
bremsstrahlung cooling can radiate energy from the cloud at 
a rate ~ r i n e n t K coo \{Z 1 T). Cooling and CR heating are thus 
approximately balancing each other for systems with a radius 

. _ JyuaPcr _ f s X C R VAkT 

Acrit- 7 7 , U-V 

n e n,A coo ] (A, n e A coo i 

where \i, =n,/n = (3Xh+ 1)/(5Xh+3) ~ 0.48 assuming primor- 
dial element composition. 

In Fig.@]we show A cr i t as a function of radius for the parame- 
ters of M87 used throu ghou t, i.e., i>a = 130 km s -1 , A"cr = 0.31, 
A C ooi(7bi-Zo) of Section [3~2l and the observed electron temper- 
ature and density profiles of Equations (f2) and (01. The CR 
heating rate without azimuthal fluctuations (/ sup = 1) yields a 
critical wavelength A cl -; t > r at all radii, which implies that the 
unstable wavelengths cannot be supported by the cluster core. 

Since CRs stream anisotropically along magnetic fields, this 
may suppress the heating flux into a given volume element de- 
pending on the magnetic connectivity to the surroundings. In 
the case of an isotropically tangled small-scale magnetic field, 
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Figure 4. Radial dependence of the critical instability length scale A cl ; t that is 
obtained by balancing the CR heating rate "Hen with the radiative cooling rate 
C ra d- On scales A < A cr j t , CR heating dominates over cooling and wipes out 
any temperature inhomogeneities. Gas parcels on scales A > A cl -j t can become 
thermally unstable unless A > r, for which the unstable wavelengths cannot 
be supported by the system. To compute the radiative cooling, we vary the 
abundance of the gas between Z = O.IZq (blue) and Z = 1.3Zq (red). We 
show one case without azimuthal fluctuations in the CR heating rate (dashed) 
and one with magnetic suppression by a factor / sup = 0.3 (solid). 

CR heating along a given direction is suppressed by a factor 
/sup ~ 0.3. As a result, the minimum size of an unstable region 
at 1 kpc is A cl it = 5-8 kpc, depending on metallicity. However, 
such a large region is unlikely to be exclusively heated from 
one direction, which makes it challenging for cooling multi- 
phase gas to form unless a given regions is sufficiently magnet- 
ically insulated with a suppression factor / sup <C 1 . Since we 
have demonstrated in Fig.|2]that CR heating globally balances 
cooling, a magnetic suppression of heating into one region im- 
plies an overheating in another. Because the gas is thermally 
unstable for temperatures in the range of (1 - 3) x 10 7 K, the ad- 
ditional energy would be radiated away on the radiative cooling 
time. 

3.5. The Emerging Physical Model of a " 'Cooling Flow" 

The previous considerations enable us to construct a compre- 
hensive "cooling flow" model for M87, which may serve as a 
blueprint for other CC clusters. We distinguish three different 
cases, each of which may be realized by a different magnetic 
suppression factor. However, quantifying their relative abun- 
dance requires appropriate numerical simulations. 

Case A. In the absence of strong magnetic suppression of the 
CR heat flux, a gas parcel is thermally stabilized by CR heating 
because the system cannot support the unstable wavelengths, 
which are larger than the cluster radius. 

Case B. For strong magnetic suppression of the CR heat flux, 
gas at temperatures T = (1-3) x 10 7 K is thermally unsta- 
ble, starts to cool radiatively, becomes denser than the ambi- 
ent gas and begins to sink toward the center. Since CRs are 
largely confined to stay on individual field lines, which are 
flux frozen into the contracting gas, the CR population is adi- 
abatically compressed. This implies an increasing CR heating 
rate Hcr oc /9 7+1//3 (Equation [Toll, which however increases at 
a slowe r rat e than the cooling rate for T = (1— 3) X 10 7 K (see 
Section [3~3l l. The CRs transfer energy to the gas on a timescale 
comparable to the Alfven crossing time over the CR pressure 
gradient length, ta = SI/va ~ 75 Myr for SI = 10 kpc. 
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If the gas parcel remains magnetically insulated, it cannot 
be supplied with a fresh influx of CRs, which causes the CR 
heating to lag the radiative cooling, implying eventually the 
formation of multi-phase gas and potentially stars, provided 
favourable conditions for their formation are met. Such a mag- 
netic suppression could be realized through processes that dy- 
namically form coherently aligned magnetic sheaths, which are 
wrapped around the gas parcel as a result of either the action 
of magnetic draping (lLvutikovll2006t iRuszkowski et alj|2007t 
iDursi & Pfrommerll2008l). the heat-flux driven b uoyancy insta- 
bility (IOuataertll2008HParrish & Ouataertl2008l) . or through ra- 
dial collapse, which amplifies the azimuthal components of the 
magnetic field, in combination with coherent motions of the 
cooling gas parcel so that the azimuthal magnetic field becomes 
causally unconnected to its surroundings. 

The comparably large critical wavelength of the instability, 
A CI it, at large radii (even for substantial magnetic suppression) 
may make it difficult for thermal instability to condense out 
large gas clouds (which requires a stable magnetic insulation 
on the correspondingly long timescale ta). Instead, it is more 
likely that cooling (multi-phase) clouds form in the central re- 
gions at radii 3 kpc because of the small A cr if For this case, 
the presence of small pockets of multi -phase gas at larger radii 
will then necessarily imply that they would have to be uplifted 
by eddies driven by the rising radio cocoons. 

Case C. For moderate magnetic suppression of a given gas 
parcel atT = (l-3)x 10 7 K, CR heating will be able to approx- 
imately match radiative cooling. That means we can consider 
the gas parcel to be in global thermal equilibrium but subject to 
local thermal instability. If this global thermal equilibrium can 
be maintained during the cooling of a necessarily large region 
of wavelength A ~ r, and if the perturbations are in pressure 
equilibrium with the surroundings, then the results of our sta- 
bility analysis applies. In particular, the gas collapse is halted 
at around 1 keV by the first "island of stability" (situated at 
Tfloor = (0.7- 1) x 10 7 K). Hence for case C, we predict a tem- 
perature "floor" at the cluster center with typical values dis- 
tributed within Tfloor. The insensitivity of the shape of the first 
"island of stability" on the various parameters entering the in- 
stability criterion suggests that we found a universal mecha- 
nism that sets a lower temperature floor in CC clusters (as long 
as the CR pressure is sufficient to provide global thermal sta- 
bility). 

This prediction of a generic lower temperature floor at Tf\ om 
is in agreement with high-resolution x-ray obs ervations of 
Virgo (M atsushita e"t~aTl l2002; Mill ion et al.ll2010l) . suggesting 
a viable solution for the puzzling absence of a large amount of 
lower-temperature gas in Virgo that is expected for an unim- 
peded cooling flow. However, gas in the high-density tail of 
density fluctuations can cool below T ~ 0.7 x 10 7 K, and hence 
leave this "island of stability". These dense fluctuations again 
become subject to therma l instability, spurring some of the ob- 
servable multi-phase gas (IWerner et al.ll20l61 120121) . The col- 
lapse could again be halted at the second and third "island of 
stability" just above 10 5 K and 10 6 K. However, since cooling 
clouds at these temperatures are associated with higher densi- 
ties and smaller length scales, the amount of CR heating now 
depends critically on the magnetic connectivity of the cloud 
with the surroundings since that determines whether there will 
be a sufficient influx of CRs responsible for sustaining the heat- 
ing rate. 

We emphasize that the success of our model came about 
without any tunable parameters and only used input from non- 



thermal radio and gamma-ray observations. The model implies 
that gas in the core region of Virgo is in global equilibrium 
but susceptible to local thermal instability. Hence this provides 
a physical foundation of the phenomenological prescription 
used in numerical simulations that study the interplay of ther- 
mal instability and feedback regulation dMcCourt et al.ll20l2t 
ISharma et al.ll2012T) . There are also notable differences since 
these prescriptions did not have the modulating effect of our 
"islands of stability". We may speculate that our model may 
then inherit some of their successes including quantitative con- 
sistency with data on multi-phase gas and star formation rates 
but holds the promise to cure the problem of the absence of a 
(physical) temperature floor in these numerical experiments. 

4. OBSERVATIONAL PREDICTIONS 

This model makes several predictions for gamma-ray, radio, 
and Sunyaev-ZeFdovich observations, which we will discuss in 
turn. Some of the proposed observations may first be possible 
in M87, owing to its immediate vicinity. To address the uni- 
versality of the this heating mechanism, these observations are 
certainly worthwhile to be pursued in other nearby CC cluster 
that show strong signs of AGN feedback. 

4.1. Gamma Rays 

If the gamma-ray emission of the low-state of M87 results 
from hadronic CRp-p interactions then the emission should be 
extended with a characteri stic angular profile worked out by 
iPfrommer & Enfilinl (120031) for their isobaric CR model. The 
gamma-ray emission should not show any sign of time variabil- 
ity and should be characterized by the distinct spectral signa- 
ture of decaying pions around 70 Me V (in the differential spec- 
trum). The triad of these observations will prove the existence 
of CRps mixed into the ICM and measure their pressure con- 
tribution. Spatially resolving the gamma-ray emission helps in 
estimating the large-scale CR pressure gradient and its contri- 
bution to the CR-heating rate. However, this does not shed light 
on the contribution from small-scale pressure fluctuations. The 
comparably coarse angular resolution of Fermi (68% contain- 
ment angle of 12' above 10 GeV) precludes such a measure- 
ment in M87. In contrast, the superior resolution of current 
and future Cerenkov telescopes (radius of a circle enclosing 
80% of the light is 1.4' with H.E.S.S. II above TeV energies 
and expected to be five times better with CTA) is a promis- 
ing route toward measuring the pressure gradient of the high- 
energy CRps, which however carry only a small amount of the 
total pressure that is dominated by the GeV-energy regime. 

We note that a critical prediction of this model is the simi- 
larity of the spectral indexes of the gamma-ray emission above 
GeV energies (which corresponds to that of the CR population) 
and the injection index of the CRe population, as probed by 
low-frequency radio observations of the core regions of Virgo 
A. Similar observations of nearby AGNs in CC clusters are 
necessary to show the universality of the presented mechanism. 

4.2. Radio Emission 

Faraday rotation measurements of background sources or 
polarized regions of Virgo A (central region, halo, or rising 
cocoons) will provide not only a reliable estimate of the spa- 
tial distribution of field strengths, but also about the magnetic 
morphology. This is important for predicting the heating rates 
(that scale with the Alfven velocity) and the microphysics of 
CR Alfven wave heating, which depends on in situ alignments 
of magnetic fields with the CR pressure gradients. In partic- 
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ular, CR release time scales from the radio cocoons may de- 
pend on the pre-existing morphology of the magnetic fields that 
dynamically drape around a super-alfve nically rising cocoon 
dLvutikovll2006tlRuszkowski et alj|2007l) . This magnetic drap- 
ing effect stabil izes the inter faces against the shear that creates 
it in first place (Dursi 2007). However, magnetic draping can 
only occur if the ratio of the integral length scale of the mag- 
netic power spectrum Xb and the curvature radius R of the radio 
cocoon at the stagnatio n line obeys the relation (see Supple- 
mentary Information of lPfrommer & Dursill2010l) 

^>M^-(m^)\o.2SM-'(1)-"\ (,4> 

where M and A4a are the sonic and Alfvenic Mach numbers 
of the rising cocoons, 7 = 5/3, and f3 is the plasma beta factor. 
The violation of this relation could cause Kelvin-Helmholtz 
and Rayleigh-Taylor instabilities to shred the cocoons and to 
release CRs into the ICM. Higher-order Faraday rotation mea- 
sure statistics are sensitive to the statistics of the projected 
magnetic field distribution. Complementary to this are mag- 
netic draping studies at spiral galaxies within this central re- 
gion (iPfrommer & Dursil 120101) that provide a means of mea- 
suring the projected in situ geometry of the magnetic field at the 
galaxies' position. Additionally, those could be of paramount 
importance to provide indications for a spatially preferred di- 
rection of magnetic fields which could hint at the presence of a 
magnetic buoyancy instability. 

The hadronic CRp-p reaction also implies the injection of 
secondary CRes that contribute to the radio synchrotron emis- 
sion with a spectral index that reflects the CRp proton distri- 
bution (i.e., ay = (a- l)/2 = .63 in M87). Using the formal- 
ism of lPfrommer et al.l (120081) . we calculate the hadronically- 
induced radio halo emission within a spherical region of radius 
34 kpc (adopting the same parameters as in the rest of the pa- 
per, namely Xqr = 0.31, a = 2.26 , q = 0.58 and the model for 
the magnetic field of Section ITU At 1.4 GHz (10.55 GHz), 
this emission component is expected to have a total flux den- 
sity of 21.1 Jy (2.15 Jy), which falls short of the observed 
ha lo emission by a factor of 6.6 (10) (see Figures 10 and 11 
of lde Gasperin et al.ll20"l2l) . 

However, the spatial profile of the expected hadronically- 
induced emission is very different from the observed radio 
halo emission, which is centrally more concentrated. If we ex- 
clude a cylinder centered on the central radio cocoon of radius 
2.5 kpc (corresponding to 0.51'), we obtain a flux density at 
1.4 GHz (10.55 GHz) for the external hadronic halo of 19.8 Jy 
(2.0 Jy), which compares favorably to the observed emission 
of 69 .1 Jy (3.4 Jy) that subten ds the same solid angle (Figure 
1 1 of lde Gasperin et aill2012l) . In particular, the hadronically- 
induced emission is expected to make up ~ 60% of the ob- 
served (external) halo emission at 10.55 GHz. Clearly, the pos- 
sibility of such an emission component should be accounted 
for when drawing conclusions about the CRe pitch angle dis- 
tribution in spectral agei ng models of synchrotro n and inverse 
Compton emitting CRes dde Gasperin et al.l2012l) . especially if 
the high-frequency data is sparse so that individual data points 
impart high statistical weights in the fits. 

If CR electrons and protons have a very different distribu- 
tion within the radio halo (which could be achieved through en- 
hanced radiative losses of CRes in regions of strong magnetic 
fields), there may be regions in which hadronically-induced 
emission prevails over the primary emission. Those regions 



should then show a radio spectral flattening in comparison to 
the adjacent primary emission-dominated regions. Because of 
the flatter slope of the hadronically-induced emission, such a 
measurement may be even more promising at high radio fre- 
quencies v >, 10 GHz. 

4.3. Sunyaev-Zel'dovich observations 

High-resolution Sunyaev-Zel'dovich studies can reveal 
the composition of the plasma within radio cocoons 
(IPfrommer et al.l 120051) . This could unambiguously demon- 
strate the presence of CRs in bubbles which would imprint as 
high-contrast cavities in high-resolution Sunyaev-Zel'dovich 
maps at the radio cocoons' positions. Their signature is dif- 
ferent from the hypothetical alternative of a hot, but sub- 
relativistic particle population. However, the large dynamic 
range required for such a detection in Virgo A will make such 
a measurement very challenging. 

5. DISCUSSION AND CONCLUSIONS 

What is the dominant physical mechanism that quenches 
cooling flows in CC clusters? Over the last decade, many 
different processes have been proposed. We will discuss the 
pros and cons of the major contenders that are representative of 
generic classes of heating models. For the first time, we have 
proposed a comprehensive CC heating model that not only con- 
vincingly matches radio, x-ray, and gamma-ray data for M87 
but is also virtually free of tunable parameters. In particular, 
our global and local stability analyses of the CR heating model 
reveal generic aspects that may apply universally to CC clus- 
ters, provided the CR release timescales from the radio cocoons 
are similar to that in M87. 

5.1. Overview of Solutions to the "Cooling Flow Problem" 

Thermal conduction of heat from the hot outer cluster re- 
gions to the cooling cores has been proposed to solve the 
prob lem (Nar avan & Medvedevl |200H IZakamska & Naravanl 
120031) . However, there are a number of drawbacks of this 
proposal. (1) The strong temperature dependence of ther- 
mal conductivity (k oc T 5 / 2 ) makes this a very inefficient pro- 
cess for groups and small clusters. (2) Thermal conduction 
ceases to be important if gas cools to low temperatures and, 
as such, requires in some clusters a conductivity above the 
clas sical Spitzer value as well as a high degree of fine tun- 
ing (IVoigt & Fabianll200l IConrov & Ostrikei1l2008l). (3) The 
heat-flux-driven buoyancy instability (Ouataert 2008) causes 
the magnetic field lines to align with the gravitational equipo- 
tential su rfaces of the cluster in the non-linear stage of the in- 
stability dParrish & Ouataertll2008l) . At first sight, this should 
suppress the inward transport of heat by more than an order 
of magnitude and should reinforce the cooling flow problem. 
The instability, however, is easily quenched by the presence of 
exter nal turbulence at the level of ~ 1% of the t hermal pres- 
sure dParrish et al.ll2010l iRuszkowski & Ohll201QT) . leaving the 
question of how much thermal conduction is suppressed by 
the presence of magnetic fields in realistic cosmological situa- 
tions. Overall, thermal conduction could certainly aid in mit- 
igating the fast onset of cooling at intermediate cluster radii 
dVoitll201 lb but is unlikely to play a major role in the cluster 
centers. 

Over the last years, AGN feedback has become the lead ac- 
tor because of its promising overall energetics and the obser- 
vational confirmation that the nuclear activity a ppears to be 
tightly coupled to a self -regulated feedback loop (iBtrzan et al.l 
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12004) . Guided by high-resolution Chandra x-ray observations, 
which emphasize sharp density features in high-pass filtered 
maps as a result of the gas being pushed around by AGN jets 
and cocoons, many papers have focused on the dissipation of 
PdV work done by the inflating cocoon on the surroundings. 
The detailed heating mechan isms include sound waves, w eak 
shocks, and turbulence (see iMcNamara & Nulsenl 120121 and 
references therein). Since sound waves and weak shocks are 
excited in the center by the AGN, their outward-directed mo- 
mentum flux in combination with the small dissipation effi- 
ciency implies that only a small fraction of this energy will be 
used for heating the central parts. Those parts, however, have 
the lowest cooling time. Additionally, the implied heating rate 
is very intermittent in space and time with an unknown duty 
cycle. Nevertheless, these processes will certainly happen and 
(if the outburst energy is large enough) may l ocally even dom- 
inate during some times dNulsen et al.ll2005l) . As we argue in 
this paper, these processes are however unlikely to be the most 
important heating source because they do not provide a generic 
solution to the lower temperature floor observed in CCs and 
they are not favored on energetic grounds due to their small 
heating efficiency. 

Assuming a relativistic filling of the cocoons, three quarters 
of the total cocoon enthalpy is stored in internal energy. We 
provide arguments in favor of a CRp population that domi- 
nates the internal cocoon pressure. Once released, those CRp 
stream at about the Alfven velocity with respect to the plasma 
rest frame and excite Alfven waves through the "streaming in- 
stability" dKulsrud & Pearcdll969h . Damping of these waves 
transfers CR energy and momentum to the thermal gas at a 
rate that scales with t he product of the A lfven velocity and the 
CR pressure gradient dGuo & Ohll2008l) . While very plausible, 
there was the open question of whether the CRp population has 
been mixed into the ICM and whether it is abundant enough to 
provide sufficiently strong heating to offset the cooling catas- 
trophe. In this paper, we answer both questions affirmatively 
for the case of M87 on the basis of a detailed analysis of non- 
thermal observations. Moreover, in addressing the issue of lo- 
cal thermal stability in this model, we put forward an attractive 
solution to the problem of an apparent temperature floor in CC 
clusters. 



5.2. Cosmic Ray Heating as the Solution? 

The puzzling absence of fossil CRes implied by the LOFAR 
observations of M87 can be naturally explained by the ob- 
served shredding of radio cocoons during their rise in the gravi- 
tational potential, e.g., through the action of Kelvin-Helmholtz 
and Rayleigh-Taylor instabilities. This releases and mixes CR 
electrons and protons into the ambient dense ICM. We argue 
that mixing is essential in this picture so that the CRes have di- 
rect contact with the dense ambient plasma, which causes CRes 
to efficiently cool through Coulomb interactions. This implies 
a maximum CRe cooling time scale of ~ 40 Myr, which is sim- 
ilar to the age of the Virgo A radio halo. A direct consequence 
of this interpretation is that CRps are similarly mixed with the 
dense gas and must collide inelastically with gas protons. This 
produces pions that decay into a steady gamma-ray signature, 
which resembles that of the flux in the "low state" of M87 ob- 
served by Fermi and H.E.S.S. In particular, the CRp population 
(witnessed through the gamma-ray emission) shares the same 
spectral index as the injected CRe population in the core region 
of Virgo A, which we consider as strong support for a com- 



mon origin of the radio and gamma-ray emission^ Hence, the 
hadronically induced gamma-ray emission enables us to esti- 
mate the CR pressure, which has a CR-to-thermal pressure con- 
tribution of Xqr = 0.31 (allowing for a conservative Coulomb 
cooling timescale of order the radio halo age that modifies the 
low-energy part of the CRp distribution). 

We show that streaming CRs heat the surrounding thermal 
plasma at a rate that balances that of radiative cooling on aver- 
age at each radius. Hence, the thermal pressure profile is a con- 
sequence of the interplay of CR heating and radiative cooling 
and should match that of the streaming CRs (hence justifying 
in retrospect our assumption of a constant CR-to-thermal pres- 
sure in our analysis). We argue that this similar radial scaling 
of the heating and cooling rate profiles is not a coincidence in 
M87 but rather generically expected in models of CR Alfven 
wave heating. 

However, the global thermal equilibrium is locally unstable 
if the gas temperature drops below kT ~ 3 keV (for solar metal- 
licity). The fate of a gas parcel depends on its magnetic con- 
nectivity to the surroundings, which may obey a yet unknown 
distribution function. Nevertheless, we can distinguish three 
different cases. 

(1) For magnetically well-connected patches without appre- 
ciable fluctuations in the CR heating rate, temperature inhomo- 
geneities are smoothed out on scales smaller than the cluster- 
centric radius. As a result, the system is stabilized because it 
cannot support the unstable wavelengths. 

(2) Contrarily, for strong magnetic suppression of the CR 
heat flux, the thermally unstable gas starts to cool radiatively, 
becomes denser than the ambient gas, and sinks toward the cen- 
ter. After an initial phase when the magnetically trapped CRs 
within the parcel deposit their energy to the gas, CR heating 
lags cooling. This eventually causes the formation of multi- 
phase gas and potentially stars. At large radii, the critical wave- 
length A cr it above which thermal instability can operate is com- 
parably large, thereby making it more likely for multi-phase 
gas to form further inward where A cr j t is smaller. Thus, small 
clumps of multi-phase gas at larger radii would have to be up- 
lifted by eddies driven by the rising radio cocoons. 

(3) For the perhaps most probable case of moderate mag- 
netic suppression, CR heating will be able to approximately 
match radiative cooling. If this global thermal equilibrium can 
be maintained during the Lagrangian evolution of the cooling 
gas parcel that is in pressure equilibrium with its surroundings, 
then the gas collapse is halted at around 1 keV by an "island of 
stability". There, the CR heating rate has a steeper dependence 
on temperature than the radiative cooling rate. Increasing (de- 
creasing) the temperature of a perturbed parcel of gas implies 
a stronger cooling (heating) that brings the parcel back to the 
background temperature. This is consistent with the observed 
temperature floor of 0.7- 1 keV seen in the x-ray observations 
of the center of Virgo. Since this result is insensitive to signif- 
icant variations of the CR pressure support, this may suggest 
a viable solution to the "cooling flow problem" (as long as the 
CR heating rate is sufficient to provide global thermal stabil- 
ity). 

However, the high-density tail of density fluctuations that is 
associated with temperatures kT <, 0.7 keV becomes again sub- 
ject to thermal instability, thereby providing a plausible mech- 

6 As|de Gasperin et al. (2012), we interpret the steeper CRe spectra in the 
halo of Virgo A as a result of a superposition of CRe spectra with different 
radiative cooling breaks that they attained during their transport from the core 
region, which is an effect that should leave the CRp spectra invariant. 
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anism for the formation of some of the observed cooling multi- 
phase medium in and around M87. We note that the success 
of the CR heating model came about without any tunable pa- 
rameters and only used input from radio and gamma-ray obser- 
vations. We present several testable predictions of this model 
that include gamma-ray, radio, and Sunyaev-Zel'dovich obser- 
vations. 

We emphasize the importance of non-thermal observations 
to establish this physical picture of AGN feedback and to pin- 
point the details of this heating mechanism in M87. While 
tempting, we caution against directly extrapolating our result 
to other CC clusters. To faithfully do so, we would need simi- 
lar non-thermal observations to establish that CRs are injected 
and mixed into the ambient ICM and to estimate the CR pres- 
sure contribution that determines the CR heating rate. In par- 
ticular, CR release time scales from the radio cocoons are cru- 
cial and may depend on the pre-existing morphology of the 
magnetic fields that dynamically drape around the rising co- 



coons dRuszkowski et al.1 120071) . Those draped fields provide 
stability in the direction of the magnetic field against the shear 
that creates it if the magnetic coherence scale is large enough 
for draping to be established in the first place. Future Fara- 
day rotation measure stud ies and galaxy draping observations 
dPfrommer & D ursi 2010) may help to clarify this point. 

We hope that this work motivates fruitful observational and 
theoretical efforts toward consolidating the presented picture 
and to establish non-thermal observations of the underlying 
high-energy processes as a driver for the understanding of the 
evolution of cluster cores and perhaps the modulation of star 
formation in elliptical galaxies. 

I thank Ewald Puchwein and Volker Springel for critical 
reading of the manuscript and helpful suggestions. I thank 
Heino Falcke and Frank Rieger for discussions and gratefully 
acknowledge financial support of the Klaus Tschira Founda- 
tion. 



APPENDIX 

A. THERMAL STABILITY IN THE PRESENCE OF A RELATIVISTIC FLUID 

A. 1. Balancing Bremsstrahlung and Line Cooling by Generic Heating Mechanisms 

Consider a fluid subject to external heating and cooling. The generalized Field criterion (Field 1965; Balbus 1986) for thermal 
instability in the presence of isobaric perturbations is given by 

d(C/T) 



dT 



<0. 



(Al) 



Here £ is the net loss function in units of ergs s" 1 g" 1 defined such that pC = C — H, where C and H are the cooling and heating 
rates per unit volume. C is generally a function of gas mass density p, tempe rature T, and possibly spa ce and time. We study the 
instability in the presence of a relativistic fluid and generalize the derivation by lConroy & Ostrikerl (120081) to also include cooling by 
line emission in addition to thermal bremsstrahlung, which will turn out to be important for putting the complete picture in place (in 
particular for cooling flows in clusters). 
A gas in thermal equilibrium (pC = 0) obeys the following identity, 
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dT 
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(A2) 



The first and third derivatives on the right-hand side are governed by the heating and cooling processes under consideration, while the 
second depends on the various sources of pressure support. We consider pressure support provided by thermal gas and a relativistic 
fluid (dominated by CRps): 

P = P th + PcR = K lP T+P C R(p), 



where K\ is a constant and Pcr is only a function of density. We use the following identity 
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(A4) 



(A5) 



where we defined Xqr = Pc R /Pth- 

We now assume that the heating has a power-law dependence on p and T and account for cooling by thermal bremsstrahlung and 
line emission to get 

P c = c rdd -n = p 2 [MTM+Aj^Zft-roffT- 1 /™, (A6) 

where A h (T,Z) = A Z 2 T^ 2 and A,(T,Z) are the cooling functions due to bremsstrahlung and line emission, respectively. Both 
depend on metallicity and temperature. (For brevity, we omit the arguments in the followin g.) Ao and To are constants, and we 
define the identical power-law indices S and j3 for the heating as in lConrov & Ostrikerl (12008b to ease comparison with their results. 
Using the relation for thermal equilibrium (pC = 0), we have 



d(pC) 



dT 



= A b p 2 T~ 
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and 



d(pQ 



dp 



1 + T b 



Combining those with Equations iA2i and ( |A5t , we finally obtain 



d(pQ 



dT 



A h pT~ 
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+ --S 
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To simplify the instability criterion, we define the logarithmic density slope of the CR pressure 7, the logarithmic temperature slope 
of the line cooling function \, and the ratio of line-to-bremsstrahlung cooling, 7?: 



<91nP, 



7: 



CR 



<91np 



X = X(T,Z) 



3 In A, 



dlnT 



A; 

and j] = j](T,Z)= — . 



Thus, the gas is thermally unstable if the discriminant T> is negative: 



V =l-5+ij 



(2-/3X1 + 7?) 



<0. 



(A10) 



(All) 



1 +X CR 7 

For 77 = we recover the result of lConrov & Ostrikerl (120081 ) if we identify Xqr = 1 /a' in their notation. Clearly, the presence of line 
cooling (?/ > 0) can substantially modify the instability criterion, depending on the temperature slope \- 

A.2. Thermal Stability of a Radiatively Cooling Gas Heated by Cosmic Ray-Excited Alfven Waves 

Now, we specify this instability criterion to our case at hand of CR Alfven wave heating. On short timescales that we consider 
for our stability analysis, the CR population can be considered adiabatically since it does not s uffer non-adiabatic losses such as 
hadronic interactions. In particular, CR streaming is an adiabatic process for the CR population (En filin et al.ll2.Ql ll) . implying that 
the CR pressure only scales with the gas density as Pcr oc p 1 yielding 

Hcr = T p^T-^ 2+s = -v a VP C r cx pi +l '\ (A12) 

Here we identified the characteristic size of our unstable clump with an exponential pressure scale height h so that the local gradient 
length scales as 51 ~ h oc p~ 1//3 and we assume va = const, as we do in the reminder of this work. In practice, va may inherit 
a weak density scaling depending on the topology of the magnetic field t hat is compressed or expanded. We defer dynamical 
simulations of this interacting CR-MHD system to future work. Equation ( IA1 2b defines the power-law indices (3 = 7 + 1 /3 and 
5=1/2 for the heating rate. We emphasize that this local adiabatic behavior of Pqr may not hold for the equilibrium distribution 
that is subject to gain and loss processes, which modify the "CR entropy" (i.e., its probability density of the microstates). In fact, 
we argued in the main part of the paper that CR Alfven wave heating establishes an equilibrium distribution that is characterized by 
Xqr = Pcr/Pi\i — const, if the dominant heating source is provided by CRs. 

To estimate the CR adiabatic exponent, we remind ourselves of the definition of the CR pressure of a CR power-law population 
defined in Equation (fTJ, which is given by 

„2 
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where v/c = p/ y 1 + p 2 is the dimensionless velocity of the CR particle and B x (a , b) denotes the incomplete beta function. The local 
adiabatic exponent of the CR population is obtained by taking the logarithmic derivative of the CR pressure at constant CR entropy 
Scr 
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1.37 



(A14) 



where we used typical values for M87 of a = 2.26 and q = 0.58 as derived in Section I2T21 

The ratio of line-to-brem sstrahlung cooling at 1 keV for a gas with solar abundance is approximately 77 ~ 3.6 (see Fig. 8 of 
iSutherland & Dopita|[T99l . assuming collisional ionization equilibrium. We can now rewrite Equation (IA1 11 to get a constraint for 
the temperature slope of the cooling function x f° r which the gas can become thermally unstable, 



1 (2-/3X1+7?) 
X 2t? (1+XcR^r, 



: 0.13, for 1 keV and Z = Z 



(A15) 



Hence, if the cooling function is sufficiently flat or even inverted (as it can be the case in the line cooling regime) we have thermal 
instability. In the bremsstrahlung regime, on the contrary, the gas is in thermal equilibrium, since the criterion for instability 



6> 



Xcr7-1 + /3 



0.8 



I+XCRJ 

cannot be fulfilled for CR Alfven wave heating that does not depend on temperature (5=1 /2). 



(A16) 
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